clear; clc; close all;

u   = 0.1;
k   = 1.8;
a   = -1;
b   = 15.2;   % 此时出现 period-7
x0  = 0.1;
q0  = 0.1;

N        = 1000000;   % 总迭代步数
Ndiscard = 1000;    % 舍弃前 500 步
Nkeep    = N - Ndiscard;  % 用于分析和绘图的步数

%% 2. 多步迭代, 返回 (xArr, qArr)
[xArr, qArr] = iterate_nhmm(u, k, a, b, x0, q0, N);

% 只取后 Nkeep 步
xPlot = xArr(end-Nkeep+1 : end);
qPlot = qArr(end-Nkeep+1 : end);

%% 3. 检测周期 7
period = 7;

% 为了安全, 在最后 14~21 步左右检测(或更多), 这里举例只检测最后 14 步
% indexEnd = length(xPlot);   % 也就是 Nkeep
% idx1 = indexEnd - period*2 + 1; % 比如最后 14 步的起始索引
% idx2 = idx1 + period - 1;   % idx1 ~ idx2 为后 7 步
%
% 这里仅做一个简单检查，如果相差很小就认为出现周期 7

idxStart = length(xPlot) - 14 + 1; % 最后14步的起始
segment1_x = xPlot(idxStart : idxStart+period-1);
segment1_q = qPlot(idxStart : idxStart+period-1);

segment2_x = xPlot(idxStart+period : idxStart+2*period-1);
segment2_q = qPlot(idxStart+period : idxStart+2*period-1);

% 计算两段的差异
diff_x = abs(segment1_x - segment2_x);
diff_q = abs(segment1_q - segment2_q);
maxDiff = max([diff_x; diff_q]);

disp(['两组点坐标最大差 = ', num2str(maxDiff)]);
if maxDiff < 1e-3
    disp('=> 检测到稳定的 period-7 吸引子');
    % 将那 7 个周期点单独标注(后 7 步)
    plot(segment2_x, segment2_q, 'ro','MarkerSize',6,'LineWidth',1.2);
    xlabel('x'); ylabel('q');
    title(['b = ', num2str(b), ' (Period-7 xiyinzi)']);
else
    disp('=> 暂未检测到严格 period-7');
end
